function [beta export] = estimate_griliches(M)

  global omega_coeffs seomega_coeffs;
  omega_coeffs = []; seomega_coeffs = [];
  
  data.firmid = M(:,1);
  data.year = M(:,2);
  data.lnVA = M(:,3);     % Value added
  data.lnY = M(:,11);     % Output
  data.lnqempl = M(:,4);  % Quality adjusted labor from Mincer
  data.Llnqempl = M(:,12);% lag    
  data.lnL = M(:,10);     % Employment
  data.L = exp(data.lnL);
  data.LlnL = M(:,13);    
  data.lnInp = M(:,5);    % Intermediate useage
  data.LlnInp = M(:,14);  
  data.Inv = M(:,6);      % Investment
  data.lnK = M(:,7);
  data.LlnK = M(:,15);
  data.Ex = M(:,8);       % Export dummy
  data.nace = M(:,9);     % Nace 2 digit code

  % Worker variables
  firmvars = 16;
  data.sex = M(:,firmvars+4);               % Gender shares
  data.ftenure = M(:,firmvars+5:firmvars+7);% Tenure shares
  data.edu = M(:,firmvars+8:firmvars+11);   % Education shares
  data.exp = M(:,firmvars+12:firmvars+15);  % Experience shares
  data.lnw = M(:,firmvars+16);              % Average wage
    
  % Which labor variable do we want to use
  data.lnlabor = data.lnL;
  data.Llnlabor = data.LlnL;
  
  % Which outcome variable to use
  %data.y = data.lnY;
  data.y = data.lnVA;
  
  data.Inv(data.Inv<=0) = NaN;

  data.n = length(data.y);
  disp('Number of firm-year obs'); 
  disp(data.n);

  % Simple OLS (used as starting values)
  X = [ones(data.n,1) data.lnlabor data.lnK];                  % Regressors
  b_ols = inv((X'*X))*X'*data.y;                       % Coeff. estimates
  se = sqrt(diag((data.y-X*b_ols)'*(data.y-X*b_ols)/(data.n-size(X,2))*inv(X'*X)));      % Standard error
  clear X;

  % ----------------------------------------------------
  % STAGE I
  % ----------------------------------------------------
  clear lambda b;
  
  Dyear = dummyvar(data.year-1995);
  data.Dyear = Dyear(:,2:end);

  K = 2 + 12 + 39;      % Number of coeffs
  startv = [.3 0 .3 .5 0 .3 .7 .8 0 .1 .2 .2]';
  beta0 = [b_ols(1); b_ols(2); startv; zeros(K-2-12,1)];
  
  global func_Xb func_iq;
  func_Xb = []; func_iq = [];
  gg = @(beta) NLS_stage1(beta,data);

  lb = ones(K,1)*-5;
  lb(2,1) = 0;
  lb(3:14,1) = -.9;
  ub = ones(K,1)*5;

  options = optimset('Display','iter','TolFun',1e-10,'TolX',1e-10,'MaxFunEvals',5e7,'MaxIter',10000);
  [bgriliches,resnorm,resid,exitflag,output,lambda,J] = lsqnonlin(gg,beta0,lb,ub,options); 

  betaL = bgriliches(2);
  b_work = bgriliches(3:14);           % Labor characteristics coeffs  
  
  % Standard error 2 
  v = data.n-numel(beta0); 
  mse = resnorm/v; 
  cov = inv(J'*J)*mse; 
  se_gril1 = sqrt(diag(cov)); 
  
  lambda = func_Xb - betaL*log(data.L.*func_iq);
  
  % ----------------------------------------------------
  % STAGE II
  % ----------------------------------------------------
  clear beta
  gmmf = @(beta) Ackerberg_gmm(beta,lambda,data);

  b0 = b_ols(3);  
  options = optimset('Display','iter','TolFun',1e-10,'TolX',1e-10,'MaxFunEvals',5e7,'MaxIter',10000); 
  [betaK,fval,exitflag,output,grad,hessian] = fminunc(gmmf,b0,options);        
  
  gmmf(betaK);
  beta = [betaL betaK b_work' omega_coeffs];
  
  % Export output
  pr = data.y - beta(1)*log(data.L.*func_iq) - beta(2)*data.lnK;
  export = [data.firmid data.year data.nace pr func_iq];
  
end
